################################################################################################
#### 运行配置文件
################################################################################################

source /public/home/xxf2019/20231121_singleMuti/config/config.sh
conda activate archr_v2

################################################################################################
#### 下机fastq文件处理
################################################################################################
## 原始的fastq文件
cp -rf /public/home/yl2019/data_result/NOA_raw_data ${work_dir}/NOA_raw_data

## 构造library文件
mkdir -p ${output_dir}/cellranger_arc

for sample in `ls ${work_dir}/NOA_raw_data | awk -F'_' '{print $2"_"$3}' | sort -u`
do
use_sample=`echo ${sample} | sed 's/_//'`
echo "fastqs,sample,library_type" > ${output_dir}/cellranger_arc/${use_sample}_library.csv

for file in `ls ${work_dir}/NOA_raw_data | grep ${sample}`
do
atac_file_len=`echo ${file} | grep _ATAC$ | wc -l`
exp_file_len=`echo ${file} | grep _Ge$ | wc -l`
if [ ${atac_file_len} -gt 0 ]
then
echo "${work_dir}/NOA_raw_data/${file},${file},Chromatin Accessibility" >> ${output_dir}/cellranger_arc/${use_sample}_library.csv
fi
if [ ${exp_file_len} -gt 0 ]
then
echo "${work_dir}/NOA_raw_data/${file},${file},Gene Expression" >> ${output_dir}/cellranger_arc/${use_sample}_library.csv
fi
done

done

<<EOF
## csv需要满足以下格式
https://support.10xgenomics.com/single-cell-multiome-atac-gex/software/pipelines/latest/using/count
#Create a libraries CSV file
For our example, the file would look as follows:

fastqs,sample,library_type
/data/blj/singlecell/rawdata/ATAC/221352/221352A_GCRA01_N_Ge,221352A_GCRA01_N_Ge,Gene Expression
/data/blj/singlecell/rawdata/ATAC/221352/221352A_GCRA01_N_ATAC,221352A_GCRA01_N_ATAC,Chromatin Accessibility
EOF

## Run cellranger-arc count，约需要3天时间
cd ${output_dir}/cellranger_arc
for sample in `ls ${work_dir}/NOA_raw_data | awk -F'_' '{print $2"_"$3}' | sort -u`
do
echo ${sample}
use_sample=`echo ${sample} | sed 's/_//'`
mkdir -p ${output_dir}/cellranger_arc/${use_sample}
${cellranger_arc} count \
--id=${use_sample} \
--reference=${tool_path}/cellranger/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 \
--libraries=${output_dir}/cellranger_arc/${use_sample}_library.csv \
--localcores=10 \
--localmem=64
done

## RNA预处理、质控、填补并保存rdata
${Rscript_archr} ${scripts_dir}/scRNA_prepocess.R \
--input_dir ${output_dir}/cellranger_arc \
--out_path ${output_dir}/seurat_raw 

## atac的输入文件
sample_list=( testis01 testis02 testis03 )
for sample in ${sample_list[@]}
do
ln -snf ${output_dir}/cellranger_arc/${sample}/outs/atac_fragments.tsv.gz ${output_dir}/cellranger_arc/${sample}_atac_fragments.tsv.gz
ln -snf ${output_dir}/cellranger_arc/${sample}/outs/atac_fragments.tsv.gz.tbi ${output_dir}/cellranger_arc/${sample}_atac_fragments.tsv.gz.tbi
done

################################################################################################
#### 结合RNA的细胞分类，计算atac按照501bp计算peak
################################################################################################
#### 构建archr对象
${Rscript_archr} ${scripts_dir}/run_archR.callPeak.R \
--rna_data_file ${input_dir}/testis_combined.Rdata \
--input_dir ${input_dir} \
--cluster_cell_file ${config_path}/cluster_celltype.csv \
--cpu 10 \
--out_path ${output_dir}/atac_res 

################################################################################################
#### 标记并展示RNA聚类结果
################################################################################################
## 注释细胞类型
${Rscript_archr} ${scripts_dir}/annotation_celltype.R \
--input_file ${input_dir}/testis_combined.Rdata \
--cell_type_file ${output_dir}/atac_res/testis_merge_all_qc_barcode_new.txt \
--out_path ${input_dir} 

## 展示细胞亚群及marker基因
${Rscript_archr} ${scripts_dir}/rna_celltype_show.R \
--input_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--out_path ${output_dir}/rna_celltype 

################################################################################################
## 质控ATAC细胞
## RNA的异常细胞也同时去掉
################################################################################################
## 去除异常的一簇细胞，该Patchytene只在testis01患者中有比例为11%（700多个细胞），远高于别的患者的2.5%（300多个），去除以后比例正常
${Rscript_archr} ${scripts_dir}/atac_cell_check.R \
--comine_data_file ${output_dir}/atac_res/testis_combined_peak.combineRNA.Rdata \
--rna_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--out_path ${output_dir}/qc_atac
## ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata
## 20852个细胞

## 继续质控，生殖细胞和体细胞分布重聚类，去除离群的细胞
## 体细胞和生殖细胞分开
${Rscript_archr} ${scripts_dir}/atac_cell_check.v2.R \
--comine_data_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--rna_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/qc_atac_v2
## 20606个细胞 = 10833胚系 + 9773体细胞
## ${output_dir}/qc_atac_v2/all/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac_v2/all/testis_combined.annotationCellType.qc.Rdata
## ${output_dir}/qc_atac_v2/germ/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac_v2/germ/testis_combined.annotationCellType.qc.Rdata
## ${output_dir}/qc_atac_v2/somatic/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac_v2/somatic/testis_combined.annotationCellType.qc.Rdata

################################################################################################
## 去除RNA和ATAC聚类异常的细胞
################################################################################################
${Rscript_archr} ${scripts_dir}/atac_cell_check.v3.R \
--comine_data_file ${output_dir}/qc_atac_v2/all/testis_combined_peak.combineRNA.qc.Rdata \
--rna_file ${output_dir}/qc_atac_v2/all/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/qc_atac_v3
## 去除185个细胞
## 20421个细胞 = 10772胚系 + 9649体细胞

################################################################################################
## 提取最后用到的细胞的barcode
${Rscript_archr} ${scripts_dir}/get_barcode_use.R \
--comine_data_file ${output_dir}/qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/bam

## 提取bam文件
for clus in `cat ${config_path}/cluster_celltype.csv | awk -F',' '{print $1}'`
do
cat ${output_dir}/bam/barcode_use.tsv | grep -w ${clus} | awk '{print $1}' > ${output_dir}/bam/${clus}.barcode_use.tsv

bam_file=${work_dir}/input/${clus}_merge_sorted.reviewRG.bam
barcode_file=${output_dir}/bam/${clus}.barcode_use.tsv

# Save the header lines
samtools view -H ${bam_file} > ${output_dir}/bam/${clus}_merge_sorted_SAM_header
# Filter alignments using filter.txt. Use LC_ALL=C to set C locale instead of UTF-8
samtools view ${bam_file} | LC_ALL=C grep -F -f ${barcode_file} > ${output_dir}/bam/${clus}_merge_sorted_SAM_SAM_body

# Combine header and body
# Convert filtered.sam to BAM format
cat ${output_dir}/bam/${clus}_merge_sorted_SAM_header \
${output_dir}/bam/${clus}_merge_sorted_SAM_SAM_body |\
samtools view -@ 8 -b -S - > ${output_dir}/bam/${clus}_merge.reviewRG.bam

## sort and index for bam files ##
samtools sort -@ 4 -o ${output_dir}/bam/${clus}_merge_sorted.reviewRG.bam ${output_dir}/bam/${clus}_merge.reviewRG.bam 
samtools index ${output_dir}/bam/${clus}_merge_sorted.reviewRG.bam

## delete intermediate files ##
rm -rf ${output_dir}/bam/${clus}_merge_sorted_SAM_header
rm -rf ${output_dir}/bam/${clus}_merge_sorted_SAM_SAM_body
rm -rf ${output_dir}/bam/${clus}_merge.reviewRG.bam
done

################################################################################################
## 得到RNA和ATAC数据在每类细胞的平均表达及开放
################################################################################################
## 分所有的、胚系和体细胞
for type in ${type_list[@]}
do
echo $type

## call peak以及鉴定motif
${Rscript_archr} ${scripts_dir}/archr_callPeak.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac_v3/${type}

## 导出atac的矩阵,每个基因算同一类细胞的均值
${Rscript_archr} ${scripts_dir}/get_genescorematrix.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac_v3/${type}

## 导出RNA的表达矩阵,每个基因算同一类细胞的均值
${Rscript_archr} ${scripts_dir}/get_expression.R \
--rna_file ${output_dir}/qc_atac_v3/${type}/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/qc_atac_v3/${type}

## 导出peak的开放矩阵,每个基因算同一类细胞的均值
${Rscript_archr} ${scripts_dir}/get_peakscorematrix.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac_v3/${type}
done

##########################################
## 保存所有的archr对象，用archr的方式
## 用于共享很关键
for type in ${type_list[@]}
do
echo $type
${Rscript_archr} ${scripts_dir}/savaArchR_results.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac_v3/${type}
done

##########################################
## 保留一份原文件的压缩版本
# -c 表示创建新的tar包，-z 表示使用gzip算法压缩，-v 表示打印详细的输出信息，-f 表示指定输出文件名。
tar -czvf ${output_dir}/qc_atac_v3.tar.gz ${output_dir}/qc_atac_v3


################################################################################################
## 既往研究已报道的差异基因列表在我们细胞类型中的GSVA
################################################################################################
for input_file in `ls ${output_dir}/qc_atac_v3/all | grep Gene | grep tsv`
do
echo ${input_file}
o1=`echo ${input_file} | sed 's/.tsv//'`

for gene_list_file in `ls ${config_path}/reportDiffGene | grep list`
do
echo ${gene_list_file}
o2=`echo ${gene_list_file} | sed 's/.list//'`

out_name=${o1}_${o2}
${Rscript_expressionTime} ${scripts_dir}/gsva_list.R \
--input_file ${output_dir}/qc_atac_v3/all/${input_file} \
--gene_list_file ${config_path}/reportDiffGene/${gene_list_file} \
--out_name ${out_name} \
--out_path ${output_dir}/celltype_plot/gsva
done
done

################################################################################################
#### 细胞整体图谱画
################################################################################################
## 展示细胞UMAP及marker基因
## rna和atac均展示
${Rscript_archr} ${scripts_dir}/rna_celltype_show.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/celltype_plot

##########################################
## 气泡图
${Rscript_archr} ${scripts_dir}/rna_celltype_dotplot.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot

##########################################
## 差异表达基因以及差异atac基因，两种生殖细胞，一种是9类另一种是3类
## 会输出满足pct的所有基因以及满足foldchange的差异表达基因
pct=0.25
logfc=1
cpu=10
type_list=(somatic germ)
for type in ${type_list[@]}
do

########
## 差异表达(生殖细胞里面),9类
${Rscript_archr} ${scripts_dir}/diff_expression.R \
--rna_file ${output_dir}/qc_atac_v3/${type}/testis_combined.annotationCellType.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--type ${type} \
--divide 0 \
--pct ${pct} \
--logfc ${logfc} \
--cpu ${cpu} \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}

## 画图，通路富集
${Rscript_archr} ${scripts_dir}/diff_go_plot.R \
--gobp_file ${output_dir}/celltype_plot/diff_expression/${type}/one_vs_other.${type}.pct_${pct}.logfc_${logfc}.GOBP.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--pct ${pct} \
--type ${type} \
--logfc ${logfc} \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}/plot

## 画图，差异表达基因热图
${Rscript_archr} ${scripts_dir}/diff_gene_heatmap.R \
--exp_file ${output_dir}/qc_atac_v3/${type}/GeneExpression.MeanByCellType.tsv \
--diff_file ${output_dir}/celltype_plot/diff_expression/${type}/one_vs_other.${type}.pct_${pct}.tsv \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}/plot

########
## 差异atac
${Rscript_archr} ${scripts_dir}/diff_atac.R \
--atac_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--divide 0 \
--type ${type} \
--cpu ${cpu} \
--out_path ${output_dir}/celltype_plot/diff_atac/${type}

## 通路富集画图
${Rscript_archr} ${scripts_dir}/diff_go_plot.R \
--gobp_file ${output_dir}/celltype_plot/diff_atac/${type}/one_vs_other.logfc_${logfc}.GOBP.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--pct ${pct} \
--type ${type} \
--logfc ${logfc} \
--out_path ${output_dir}/celltype_plot/diff_atac/${type}/plot
done

########
## 生殖细胞分为3类
type=germ
## 表达
${Rscript_archr} ${scripts_dir}/diff_expression.R \
--rna_file ${output_dir}/qc_atac_v3/${type}/testis_combined.annotationCellType.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--type ${type} \
--divide 3 \
--pct ${pct} \
--logfc ${logfc} \
--cpu ${cpu} \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}_combine

## 画图，通路富集
${Rscript_archr} ${scripts_dir}/diff_go_plot.R \
--gobp_file ${output_dir}/celltype_plot/diff_expression/${type}_combine/one_vs_other.${type}.pct_${pct}.logfc_${logfc}.GOBP.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--pct ${pct} \
--type ${type} \
--logfc ${logfc} \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}_combine/plot

## 画图，差异表达基因热图
${Rscript_archr} ${scripts_dir}/diff_gene_heatmap.R \
--exp_file ${output_dir}/qc_atac_v3/${type}/GeneExpression.MeanByCellType.Combine.tsv \
--diff_file ${output_dir}/celltype_plot/diff_expression/${type}_combine/one_vs_other.${type}.pct_${pct}.tsv \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}_combine/plot

## atac
${Rscript_archr} ${scripts_dir}/diff_atac.R \
--atac_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--divide 3 \
--type ${type} \
--cpu ${cpu} \
--out_path ${output_dir}/celltype_plot/diff_atac/${type}_combine

## 通路富集
${Rscript_archr} ${scripts_dir}/diff_go_plot.R \
--gobp_file ${output_dir}/celltype_plot/diff_atac/${type}_combine/one_vs_other.logfc_${logfc}.GOBP.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--pct ${pct} \
--type ${type} \
--logfc ${logfc} \
--out_path ${output_dir}/celltype_plot/diff_atac/${type}_combine/plot


##########################################
## 细胞质控图
${Rscript_archr} ${scripts_dir}/cell_qc_plot.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}//qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/cell_qc

## 质控图不分样本
${Rscript_archr} ${scripts_dir}/cell_qc_plot.combine.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}//qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/cell_qc/combineSample

<<EOF
## 卞师兄质控参数
GCN2 <- subset(x = GCN,subset = nFeature_RNA>200&nFeature_RNA<4000) #10852
GCN3 <- subset(x = GCN2,subset = nCount_RNA>500&nCount_RNA<20000) #7347
GCN4 <- subset(x = GCN3,subset = nCount_ATAC>1000&nCount_ATAC<20000) #5868
GCN5 <- subset(x = GCN4,subset = nucleosome_signal < 4) #5868
GCN6 <- subset(x = GCN5,subset = TSS.enrichment > 2) #5860
GCN1 <- subset(x = GCN6,percent.mt<50) #5592

#### ATAC质控
## 目前最靠谱的质控标准是TSS富集的得分(评估ATAC-seq数据的信噪比)和唯一比对数(比对数如果不够，那么该细胞也没有分析的价值)
## https://www.archrproject.com/bookdown/per-cell-quality-control.html
## 在 ArchR 中，通过过滤细胞被识别为 TSS 富集分数大于 4 且独特核片段超过 1000 个的细胞

## 我们的数据分布
## TSS > 4 :TSS 富集评分指标背后的想法是，与其他基因组区域相比，ATAC-seq 数据在基因 TSS 区域普遍富集，这是由于与启动子结合的大型蛋白质复合物。
## nFragments below 1,000 独特核片段的数量（即未映射到线粒体 DNA）。很简单 - 具有很少可用碎片的细胞将无法提供足够的数据来做出有用的解释，因此应该被排除。
## FRIP:The fraction of reads in called peak regions, 默认不进行质控,我们样本中最低为0.148
## 质控标准TSS > 4, nFragments > 1,000

#### RNA质控
## 皮肤癌的NG文章的标准
## First, cells were removed if they had fewer than 200 genes expressed, 
## fewer than 1,000 unique sequenced reads (unique molecular identifiers) or 
## greater than 20% of counts corresponding to mitochondrial genes. 

## 我们的数据分布
## percent.mito [0,0.25]
## nCount_RNA [591,49943]:nCount_RNA is the total number of molecules detected within a cell
## nFeature_RNA [501.00,9915.00]:nFeature_RNA is the number of genes detected in each cell
## 质控标准：线粒体<0.25, 50000 > nCount_RNA > 500, 10000 > nFeature_RNA > 500
EOF

##########################################
## 展示差异的peak及其所在的motif
## 提取所有的peak及其所属类型、提取差异的peak其所属类型、提取所有存在motif的peak、提取所有motif的位置
for type in ${type_list[@]}
do
echo $type
${Rscript_archr} ${scripts_dir}/archr_markerMotif.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--divide 0 \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/diff_peak/${type}_test
done

## germ分成三大类
type=germ
${Rscript_archr} ${scripts_dir}/archr_markerMotif.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--divide 3 \
--out_path ${output_dir}/celltype_plot/diff_peak/${type}_combine

## 提取差异peak
germ_type=( "SSC_SPG" "SPT" "SPC"  )
for germ_t in ${germ_type[@]}
do
cat ${output_dir}/celltype_plot/diff_peak/${type}_combine/plotMarkerPeakHeatmap.tsv | grep -w ${germ_t} | awk '{print $1}' | tr ':-' '\t' \
> ${output_dir}/celltype_plot/diff_peak/${type}_combine/${germ_t}.peak.bed
done

<<EOF
##########################################
## 计算每个基因自身的表达及其开放程度的相关性
for type in ${type_list[@]}
do
echo $type
${Rscript_archr} ${scripts_dir}/get_exp-score_peaklink.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--cluster ${type} \
--mean_expr_file ${output_dir}/qc_atac_v3/${type}/GeneExpression.MeanByCellType.tsv \
--out_path ${output_dir}/celltype_plot/exp_score
done
EOF

################################################################################################
#### 细胞peak-gene调控关系的图
################################################################################################
##########################################
## 所有样本中peak-gene分层哪集簇、每簇对应的peak富集在哪些TF、每簇基因富集在哪些GO通路
nclust_list=`seq 11 13`
for type in ${type_list[@]}
do
echo $type

for nclust in ${nclust_list[@]}
do
echo ${nclust}
${Rscript_archr} ${scripts_dir}/archr_peak2gene.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--nclust ${nclust} \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}

## 通路富集画图
## 分前10%的peak数量的基因，以及gene其本身开放-表达显著正相关且存在至少3个peak连接的基因
${Rscript_archr} ${scripts_dir}/archr_peak2gene_go.R \
--peak_gene_file ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k${nclust}.peakNum.tsv \
--exp_score_file ${output_dir}/celltype_plot/exp_score/${type}_exp-atac_linkPeakNum.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--nclust ${nclust} \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}
done
done

## somatic的按照9的cluster分成不同peak
type=somatic
nclust=9
for clus in `cat ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k${nclust}.tsv | awk -F'\t' '{print $2}' | sort -u | grep -v kclust`
do
echo ${clus}
cat ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k${nclust}.tsv | \
awk -F'\t' '{OFS="\t"}{if($2==clus){print $10,$11,$12}}' clus=${clus} | sed 's/"//g' \
> ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k${nclust}.${clus}.bed
done

## germ的按照12的cluster分成不同peak
type=germ
nclust=12
for clus in `cat ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k${nclust}.tsv | awk -F'\t' '{print $2}' | sort -u | grep -v kclust`
do
echo ${clus}
cat ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k${nclust}.tsv | \
awk -F'\t' '{OFS="\t"}{if($2==clus){print $10,$11,$12}}' clus=${clus} | sed 's/"//g' \
> ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k${nclust}.${clus}.bed
done

## 计算所有peak的保守性，比较存在调控的peak和其它的peak的差异
${Rscript_archr} ${scripts_dir}/peak_conservation.R \
--somatic_peak_file ${output_dir}/qc_atac_v3/somatic/testis_combined_peak.qc.tsv \
--somatic_peak_peak_gene_file ${output_dir}/celltype_plot/peak2gene/somatic/peakToGeneHeatmap_LabelClust_k9.tsv \
--germ_peak_file qc_atac_v3/germ/testis_combined_peak.qc.tsv \
--germ_peak_peak_gene_file ${output_dir}/celltype_plot/peak2gene/germ/peakToGeneHeatmap_LabelClust_k12.tsv \
--out_path ${output_dir}/celltype_plot/peak2gene


<<EOF
##### 按照protein_coding 和 lncRNA分别做
## 难点在于构造表达矩阵
## 所有样本中peak-gene分层哪集簇、每簇对应的peak富集在哪些TF、每簇基因富集在哪些GO通路
type_list=(germ somatic)
gene_list=(protein_coding lncRNA protein_coding2lncRNA)
for type in ${type_list[@]}
do
for gene_type in ${gene_list[@]}
do
echo $type
echo ${gene_type}
${Rscript_archr} ${scripts_dir}/archr_peak2gene.chooseGene.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--rna_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined.annotationCellType.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--gtf_file ~/ref/GTF/20240317_gencode_v32_gene.bed \
--gff3_file ~/ref/GTF/gencode.v32.annotation.gff3 \
--gene_type ${gene_type} \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}/${gene_type}
done
done

## lcnRNA不做通路富集
nclust_list=`seq 5 20`
type_list=(germ somatic)
gene_list=(protein_coding protein_coding2lncRNA)
for type in ${type_list[@]}
do
for nclust in ${nclust_list[@]}
do
for gene_type in ${gene_list[@]}
do
echo $type
echo ${nclust}
echo ${gene_type}
## 通路富集画图
## 分钱10%的peak数量的基因，以及gene其本身开放-表达显著正相关且存在至少3个peak连接的基因
${Rscript_archr} ${scripts_dir}/archr_peak2gene_go.R \
--peak_gene_file ${output_dir}/celltype_plot/peak2gene/${type}/${gene_type}/peakToGeneHeatmap_LabelClust_k${nclust}.peakNum.tsv \
--exp_score_file ${output_dir}/celltype_plot/exp_score/${type}_exp-atac_linkPeakNum.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--nclust ${nclust} \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}/${gene_type}
done
done
done
EOF

##########################################
## 选取peak-gene中，正向相关以及不相关的peak，分别计算其上的motif和peak在每个细胞的开放程度
type=germ
peak_type_list=(positive other)
for peak_type in ${peak_type_list[@]}
do
echo ${peak_type}
${Rscript_archr} ${scripts_dir}/get_peak-gene.motif.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--rna_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined.annotationCellType.qc.Rdata \
--gene_ensg_file ~/ref/GTF/gencode.v32.gene_ensg.tsv \
--gff3_file ~/ref/GTF/gencode.v32.annotation.gff3 \
--motif_annotation_file ${output_dir}/celltype_plot/diff_peak/germ/AllPeak.containMotif.tsv \
--peak_type ${peak_type} \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}_mfuzz
done

<<EOF
## 时序聚类
## 聚成不同类
type=germ
for peak_type in ${peak_type_list[@]}
do
for nclust in `seq 5 35`
do
## peak
data_type=peak.${peak_type}
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.R \
--data_type ${data_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}_mfuzz/${data_type} \
--rsem_file ${output_dir}/celltype_plot/peak2gene/${type}_mfuzz/PeakScore.Peak-Gene.${peak_type}.rds \
--motif_file ${output_dir}/qc_atac_v3/${type}/Motif.All.rds

## motif
data_type=motif.${peak_type}
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.R \
--data_type ${data_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}_mfuzz/${data_type} \
--rsem_file ${output_dir}/celltype_plot/peak2gene/${type}_mfuzz/Motif.Peak-Gene.${peak_type}.rds \
--motif_file ${output_dir}/qc_atac_v3/${type}/Motif.All.rds
done
done
EOF

##########################################
## 提取peak，分为总的peak以及peak-gene的peak，其它的peak
## https://meme-suite.org/meme/tools/meme-chip
mkdir -p ${output_dir}/celltype_plot/peak_region
type_list=(germ somatic)
for type in ${type_list[@]}
do
echo ${type}

## 所有的peak
cat ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.qc.tsv | sed '1d' | awk -F'\t' '{OFS="\t"}{print $1,$2,$3}' | \
${bedtools} sort -i - \
> ${output_dir}/celltype_plot/peak_region/${type}.all.bed

## peak2gene的peak
cat ${output_dir}/celltype_plot/peak2gene/${type}/peakToGeneHeatmap_LabelClust_k25.tsv | \
sed '1d' | awk -F'\t' '{print $1}' | sort -u | sed 's/"//g' | tr ':-' '\t' | \
${bedtools} sort -i - \
> ${output_dir}/celltype_plot/peak_region/${type}.peak2gene.bed

## 其它的peak
${bedtools} subtract -a ${output_dir}/celltype_plot/peak_region/${type}.all.bed -b ${output_dir}/celltype_plot/peak_region/${type}.peak2gene.bed \
> ${output_dir}/celltype_plot/peak_region/${type}.other.bed
done

## 在线运行MEME的SEA

##########################################
## 计算peak在encode数据库中的富集情况
mkdir -p ${output_dir}/celltype_plot/peak_region/encode_enrich
echo -e "name1\tname2\tuniq_coverage1\tuniq_coverage2\ttotal_coverage1\ttotal_coverage2" \
> ${output_dir}/celltype_plot/peak_region/encode_enrich/region_enrich.tsv

for bedfile in `ls ${output_dir}/celltype_plot/peak_region/ | grep bed`
do
echo ${bedfile}
name1=`echo ${bedfile} | sed 's/.bed//'`
file1=${output_dir}/celltype_plot/peak_region/${bedfile}

## 文件1的覆盖区域
total_coverage1=$(awk '{sum+=$3-$2} END {print sum}' ${file1})

for encodefile in `ls ${work_dir}/encode_bed/ | grep hg38 | grep -v unmap`
do
name2=`echo ${encodefile} | sed 's/.hg38.bed//'`
file2=${work_dir}/encode_bed/${encodefile}

## 文件2的覆盖区域
total_coverage2=$(awk '{sum+=$3-$2} END {print sum}' ${file2})

# 计算第一个BED文件中的独立区域长度
uniq_coverage1=`${bedtools} subtract -a ${file1} -b ${file2} | awk '{sum+=$3-$2} END {print sum}'`
# 计算第二个BED文件中的独立区域长度
uniq_coverage2=`${bedtools} subtract -a ${file2} -b ${file1} | awk '{sum+=$3-$2} END {print sum}'`

echo -e "${name1}\t${name2}\t${uniq_coverage1}\t${uniq_coverage2}\t${total_coverage1}\t${total_coverage2}" \
>> ${output_dir}/celltype_plot/peak_region/encode_enrich/region_enrich.tsv
done
done

## 计算不同peak，其在不同encode数据的富集系数
${Rscript_archr} ${scripts_dir}/encode_region_enrich.R \
--input_file ${output_dir}/celltype_plot/peak_region/encode_enrich/region_enrich.tsv \
--out_path ${output_dir}/celltype_plot/peak_region/encode_enrich


##########################################
## 所有基因的聚类分析
## 计算motif的开放程度和tf表达的相关性，用archr自带的算法
type=germ
${Rscript_archr} ${scripts_dir}/computeCor_motif_atac-rna.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/celltype_plot/mfuzz

## 导出基因的表达矩阵，motif的开放矩阵及每类细胞的平均表达
## motif名字去除_，只保留既存在开放又存在表达的769个motif
type=germ
${Rscript_archr} ${scripts_dir}/get_expression-motif.allcell.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--rna_file ${output_dir}/qc_atac_v3/${type}/testis_combined.annotationCellType.qc.Rdata \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--out_path ${output_dir}/qc_atac_v3/${type}

<<EOF
## 趋势性聚类
## 聚成不同类
for nclust in `seq 5 35`
do
## 表达
data_type=exp
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.R \
--data_type ${data_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/mfuzz \
--rsem_file ${output_dir}/qc_atac_v3/${type}/GeneExpression.All.rds \
--motif_file ${output_dir}/qc_atac_v3/${type}/Motif.All.rds

## 所有motif的时序分析
data_type=motif
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.R \
--data_type ${data_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/mfuzz \
--rsem_file ${output_dir}/qc_atac_v3/${type}/Motif.All.rds \
--motif_file ${output_dir}/qc_atac_v3/${type}/Motif.All.rds
done

## 在不同的类中，注释TF的表达和motif的相关性，并计算显著相关的转录因子是否富集在某一类
## 不同类的基因去做gsva
cor_t_list=(0.1 0.2 0.3 0.4 0.5)
for cor_t in ${cor_t_list[@]}
do
for nclust in `seq 5 35`
do
## 表达
data_type=exp
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_positveTFenrich.R \
--rna_file ${output_dir}/qc_atac_v3/germ/GeneExpression.MeanByCellType.tsv \
--data_type ${data_type} \
--clust ${nclust} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--out_path ${output_dir}/celltype_plot/mfuzz \
--mfuzz_file ${output_dir}/celltype_plot/mfuzz/mfuzz_plot.${data_type}.${nclust}.tsv \
--cor_t ${cor_t}

## 所有motif的时序分析
data_type=motif
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_positveTFenrich.R \
--rna_file ${output_dir}/qc_atac_v3/germ/Motif.MeanByCellType.tsv \
--data_type ${data_type} \
--clust ${nclust} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--out_path ${output_dir}/celltype_plot/mfuzz \
--mfuzz_file ${output_dir}/celltype_plot/mfuzz/mfuzz_plot.${data_type}.${nclust}.tsv \
--cor_t ${cor_t}
done
done
EOF

<<EOF
#######
## 热图聚类所有基因，标记基因所属的类别，并按照类别进行富集分析
## 主要产生这种图pheatmap_exp.cluster.pdf
data_type=exp
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.R \
--data_type ${data_type} \
--rna_file ${output_dir}/qc_atac_v3/germ/GeneExpression.MeanByCellType.tsv \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--out_path ${output_dir}/celltype_plot/mfuzz
## 阳性motif富集
cor_t_list=(0.1 0.2 0.3 0.4 0.5)
for cor_t in ${cor_t_list[@]}
do
for nclust in `seq 5 35`
do
data_type=exp
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene_positveTFenrich.R \
--data_type ${data_type} \
--clust ${nclust} \
--out_path ${output_dir}/celltype_plot/mfuzz \
--cluster_file ${output_dir}/celltype_plot/mfuzz/pheatmap_${data_type}.cluster.${nclust}.gene.tsv \
--cor_t ${cor_t}
done
done
EOF

<<EOF
##########################################
## 以下分析针对motif的表达及开放
## 所有的motif表达及开放
geneset_type=all_motif
geneset_file=${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv
for nclust in `seq 5 20`
do
## 表达
data_type=exp
rsem_file=${output_dir}/qc_atac_v3/${type}/GeneExpression.All.rds
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type} \
--rsem_file ${rsem_file} \
--motif_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file}

## 所有motif的时序分析
data_type=motif
rsem_file=${output_dir}/qc_atac_v3/${type}/Motif.All.rds
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type} \
--rsem_file ${rsem_file} \
--motif_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file}
done


## 已知的motif表达
geneset_type=known_motif
geneset_file=${config_path}/Human_reported_TF2.csv
for nclust in `seq 5 20`
do
## 表达
data_type=exp
rsem_file=${output_dir}/qc_atac_v3/${type}/GeneExpression.All.rds
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type} \
--rsem_file ${rsem_file} \
--motif_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file}

## 所有motif的时序分析
data_type=motif
rsem_file=${output_dir}/qc_atac_v3/${type}/Motif.All.rds
${Rscript_expressionTime} ${scripts_dir}/time_Mfuzz_v2.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--clust ${nclust} \
--cluster_file ${config_path}/cluster_celltype.csv \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type} \
--rsem_file ${rsem_file} \
--motif_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file}
done

## pheatmap
type=germ
geneset_type=all_motif
geneset_file=${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv

## 表达
data_type=exp
rna_file=${output_dir}/qc_atac_v3/germ/GeneExpression.MeanByCellType.tsv
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--rna_file ${rna_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type}
## 所有motif的时序分析
data_type=motif
rna_file=${output_dir}/qc_atac_v3/germ/Motif.MeanByCellType.tsv
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--rna_file ${rna_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type}

## pheatmap
geneset_type=known_motif
geneset_file=${config_path}/Human_reported_TF2.csv
## 表达
data_type=exp
rna_file=${output_dir}/qc_atac_v3/germ/GeneExpression.MeanByCellType.tsv
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--rna_file ${rna_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type}
## 所有motif的时序分析
data_type=motif
rna_file=${output_dir}/qc_atac_v3/germ/Motif.MeanByCellType.tsv
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.R \
--data_type ${data_type} \
--geneset_type ${geneset_type} \
--rna_file ${rna_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/mfuzz_${geneset_type}
EOF

##################################################
## 已知重要的转录因子的motif、表达以及表达-motif的热图
geneset_type=known_motif
geneset_file=${config_path}/Human_reported_TF2.new.csv

rna_file=${output_dir}/qc_atac_v3/germ/GeneExpression.MeanByCellType.tsv
motif_file=${output_dir}/qc_atac_v3/germ/Motif.MeanByCellType.tsv

${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.combineRNA_Motif.R \
--geneset_type ${geneset_type} \
--rna_file ${rna_file} \
--motif ${motif_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/pheatmap

## 根据研究的3个level的tf
for geneset_type in `seq 1 3`
do
cat ${config_path}/Human_reported_TF2.new.csv | grep -w ${geneset_type} \
> ${output_dir}/celltype_plot/pheatmap/${geneset_type}.list
geneset_file=${output_dir}/celltype_plot/pheatmap/${geneset_type}.list
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.combineRNA_Motif.R \
--geneset_type "level"${geneset_type} \
--rna_file ${rna_file} \
--motif ${motif_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/pheatmap
done

## positve和other的peak上的motif
peak_type_list=(positive other)
for peak_type in ${peak_type_list[@]}
do
motif_file=${output_dir}/celltype_plot/peak2gene/germ_mfuzz/Motif.Peak-Gene.MeanByCellType.${peak_type}.tsv

geneset_type=known_motif
geneset_file=${config_path}/Human_reported_TF2.new.csv
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.combineRNA_Motif.R \
--geneset_type ${geneset_type}".${peak_type}" \
--rna_file ${rna_file} \
--motif ${motif_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/pheatmap
## 3个level的tf
for geneset_type in `seq 1 3`
do
cat ${config_path}/Human_reported_TF2.new.csv | grep -w ${geneset_type} \
> ${output_dir}/celltype_plot/pheatmap/${geneset_type}.list

geneset_file=${output_dir}/celltype_plot/pheatmap/${geneset_type}.list
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.combineRNA_Motif.R \
--geneset_type "level"${geneset_type}".${peak_type}" \
--rna_file ${rna_file} \
--motif ${motif_file} \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/pheatmap
done
done

################################################################################################
#### 拟时序,对鉴定的160个生殖细胞里面的TF，进行表达和motif活性的排序
################################################################################################
<<EOF
## 基于RNA
## 基于表达数据建立拟时序
${Rscript_archr} ${scripts_dir}/run_monocle.R \
--input_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/monocole

## 整合基因的atac、rna以及拟时序的结果为Rdata
${Rscript_archr} ${scripts_dir}/run_combine_rna-atac.addPseudotime.R \
--rna_data_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \ \
--rna_data_monocle_file ${output_dir}/monocole/testis.monocle.Rdata \
--atac_data_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata\
--time_diff_gene_file ${output_dir}/monocole/pseudotime_ordergene.tsv \
--gff3_file ~/ref/GTF/gencode.v32.annotation.gff3 \
--gene_ensg_file ~/ref/GTF/gencode.v32.gene_ensg.tsv \
--out_path ${output_dir}/monocole

## rna拟时序画图
## RNA本身的以及，与时序相关基因的表达和峰开放
${Rscript_archr} ${scripts_dir}/run_monocle_plot.R \
--rna_data_monocle_file ${output_dir}/monocole/testis.monocle.Rdata \
--comine_data_file ${output_dir}/monocole/testis_combined_peak.combineRNA.Rdata \
--time_diff_gene_file ${output_dir}/monocole/pseudotime_ordergene.tsv \
--out_path ${output_dir}/monocole
EOF

################################
## 基于ATAC的时序分析
## 最后使用
${Rscript_archr} ${scripts_dir}/run_trajectory.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/trajectory

################################
## 展示所有感兴趣基因集的motif
geneset_type=known_motif
geneset_file=${config_path}/Human_reported_TF2.new.csv

${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_file ${geneset_file} \
--geneset_type ${geneset_type} \
--out_path ${output_dir}/celltype_plot/trajectory/all/${geneset_type}

## 3个level的tf
for geneset_type in `seq 1 3`
do

cat ${config_path}/Human_reported_TF2.new.csv | grep -w ${geneset_type} \
> ${output_dir}/celltype_plot/trajectory/all/${geneset_type}.list
geneset_file=${output_dir}/celltype_plot/trajectory/all/${geneset_type}.list
${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_type "level"${geneset_type} \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/trajectory/all/"level"${geneset_type}
done

################################
## 提取peak-gene的peak，构建archr对象
${Rscript_archr} ${scripts_dir}/get_peak-gene.archr.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--peak_gene_file ${output_dir}/celltype_plot/peak2gene/germ/peakToGeneHeatmap_LabelClust_k25.tsv \
--out_path ${output_dir}/qc_atac_v3/germ_peak-gene

## positive里面的
geneset_type=known_motif
geneset_file=${config_path}/Human_reported_TF2.new.csv

${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_file ${geneset_file} \
--geneset_type ${geneset_type} \
--out_path ${output_dir}/celltype_plot/trajectory/positive/${geneset_type}

## 3个level的tf
for geneset_type in `seq 1 3`
do

cat ${config_path}/Human_reported_TF2.new.csv | grep -w ${geneset_type} \
> ${output_dir}/celltype_plot/trajectory/positive/${geneset_type}.list
geneset_file=${output_dir}/celltype_plot/trajectory/positive/${geneset_type}.list
${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_type "level"${geneset_type} \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/trajectory/positive/"level"${geneset_type}
done


<<EOF
## Homer只能用JASPAR的数据库，该数据库不分人类，只有脊椎动物这个大类
## motif富集分析,homer的JASPAR
mkdir -p ${output_dir}/celltype_plot/peak_region/homer
ls ${output_dir}/celltype_plot/peak_region/ | grep bed | grep peak.bed | xargs -P 10 -i sh -c '
echo {}
mkdir -p ${output_dir}/celltype_plot/peak_region/homer/{}
${findMotifsGenome} ${output_dir}/celltype_plot/peak_region/{} hg38 ${output_dir}/celltype_plot/peak_region/homer/{} -len 8,10,12
'
EOF

################################
## 计算每个基因在多少比例的细胞大于0，总的，不同的细胞中的，用于确定重要TF的界值
type=germ
${Rscript_archr} ${scripts_dir}/get_gene-expression_ratio.R \
--rna_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined.annotationCellType.qc.Rdata \
--geneset_file ${config_path}/Human_reported_TF2.addInfo.csv \
--magic_exp_file ${output_dir}/qc_atac_v3/germ/GeneExpression.MeanByCellType.magic.tsv \
--out_path ${output_dir}/qc_atac_v3/${type}/


################################
## 基于SEA对positive peak的富集结果，同时卡表达细胞的比例
geneset_type=all_motif

## 按照不同的PCT标准
for pct_type in ` cat ${config_path}/sea_exp_ratio.csv  | head -1 | tr ',' '\n' | grep pct `
do

## 判断该细胞比例阈值属于第几列
ncol=`cat ${config_path}/sea_exp_ratio.csv  | head -1 | tr ',' '\n' | grep -n -w ${pct_type} | awk -F':' '{print $1}'`
## 
cat ${config_path}/sea_exp_ratio.csv | awk -F',' '{if($ncol=="TRUE"){print$1}}' ncol=${ncol} \
> ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.list

geneset_file=${output_dir}/celltype_plot/trajectory/positive/${pct_type}.list
geneset_type=${pct_type}

${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_file ${geneset_file} \
--geneset_type ${geneset_type} \
--out_path ${output_dir}/celltype_plot/trajectory/positive/${geneset_type}
done

################################
## 基于SEA对positive peak的富集结果，同时卡表达细胞的比例，分成强相关，弱相关，不相关
geneset_type=all_motif
pct_type=pct_0.25

## 按照不同的相关系数
#cor_hh=0.7
cor_h=0.5
cor_l=0.2
cor_n=-0.2

## 强强相关
#gene_set_cor_hh=`cat ~/20231121_singleMuti/results/celltype_plot/peak2gene/germ_mfuzz/cor.motif_atac-rna.positive.tsv | \
#awk '{if($3 > cor_h)print $1}' cor_h=${cor_hh} | sed 's/"//g' | tr '\n' '|' | sed 's/|$//'`
## 强相关
gene_set_cor_h=`cat ~/20231121_singleMuti/results/celltype_plot/peak2gene/germ_mfuzz/cor.motif_atac-rna.positive.tsv | \
awk '{if($3 > cor_h)print $1}' cor_h=${cor_h} | sed 's/"//g' | tr '\n' '|' | sed 's/|$//'`
## 弱相关
gene_set_cor_l=`cat ~/20231121_singleMuti/results/celltype_plot/peak2gene/germ_mfuzz/cor.motif_atac-rna.positive.tsv | \
awk '{if($3 > cor_l)print $1}' cor_l=${cor_l} | sed 's/"//g' | tr '\n' '|' | sed 's/|$//'`
## 负相关
gene_set_cor_n=`cat ~/20231121_singleMuti/results/celltype_plot/peak2gene/germ_mfuzz/cor.motif_atac-rna.positive.tsv | \
awk '{if($3 < cor_n)print $1}' cor_n=${cor_n} | sed 's/"//g' | tr '\n' '|' | sed 's/|$//'`

## 强强相关
#cat ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.list | grep -E -w ${gene_set_cor_hh} \
#> ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.cor_${cor_hh}.list
## 强相关
cat ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.list | grep -E -w ${gene_set_cor_h}  \
> ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.cor_${cor_h}.list
## 弱相关
cat ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.list | grep -E -w ${gene_set_cor_l} | grep -v -E -w ${gene_set_cor_h} \
> ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.cor_${cor_l}.list
## 负相关
cat ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.list | grep -E -w ${gene_set_cor_n} \
> ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.cor_${cor_n}.list
## 不相关
cat ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.list | grep -v -E -w ${gene_set_cor_l} | grep -v -E -w ${gene_set_cor_n} \
> ${output_dir}/celltype_plot/trajectory/positive/${pct_type}.cor_no.list

## 画图
for geneset_file in `ls ${output_dir}/celltype_plot/trajectory/positive/ | grep cor_ | grep list `
do
echo ${geneset_file}
geneset_type=`echo ${geneset_file} | sed 's/.list//'`
geneset_file=${output_dir}/celltype_plot/trajectory/positive/${geneset_file}

${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_file ${geneset_file} \
--geneset_type ${geneset_type} \
--out_path ${output_dir}/celltype_plot/trajectory/positive/${geneset_type}
done

## 所有160个TF，画表达和motif的曲线图
geneset_file=${output_dir}/celltype_plot/trajectory/positive/pct_0.25.list
${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.smoothline.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/trajectory/positive/pct_0.25/smoothline


################################################################################################
#### 对于相关系数>0.5和大于0.2的画热图和拟时序图

## 基于peak2gene的结果
rna_file=${output_dir}/qc_atac_v3/germ_peak-gene/GeneExpression.MeanByCellType.tsv
motif_file=${output_dir}/celltype_plot/peak2gene/germ_mfuzz/Motif.Peak-Gene.MeanByCellType.positive.tsv
cor_file=${output_dir}/celltype_plot/peak2gene/germ_mfuzz/cor.motif_atac-rna.positive.tsv

#################
cor=b0.5
mkdir -p ${output_dir}/celltype_plot/pheatmap_traj_${cor}
cp -rf ${output_dir}/celltype_plot/trajectory/positive/pct_0.25.cor_0.5.list ${output_dir}/celltype_plot/pheatmap_traj_${cor}/gene.list

geneset_file=${output_dir}/celltype_plot/pheatmap_traj_${cor}/gene.list
geneset_type=positiveTF_${cor}

## 热图
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.combineRNA_Motif.R \
--geneset_type ${geneset_type} \
--rna_file ${rna_file} \
--motif ${motif_file} \
--cor_file ${cor_file} \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/pheatmap_traj_${cor}
## 拟时序
${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_file ${geneset_file} \
--geneset_type ${geneset_type} \
--out_path ${output_dir}/celltype_plot/pheatmap_traj_${cor}

#################
## 大于0.2
cor=b0.2
mkdir -p ${output_dir}/celltype_plot/pheatmap_traj_${cor}
cat ${output_dir}/celltype_plot/trajectory/positive/pct_0.25.cor_0.2.list ${output_dir}/celltype_plot/trajectory/positive/pct_0.25.cor_0.5.list \
> ${output_dir}/celltype_plot/pheatmap_traj_${cor}/gene.list

geneset_file=${output_dir}/celltype_plot/pheatmap_traj_${cor}/gene.list
geneset_type=positiveTF_${cor}

## 热图
${Rscript_expressionTime} ${scripts_dir}/pheatmap_allgene.geneSet.combineRNA_Motif.R \
--geneset_type ${geneset_type} \
--rna_file ${rna_file} \
--motif ${motif_file} \
--cor_file ${cor_file} \
--geneset_file ${geneset_file} \
--out_path ${output_dir}/celltype_plot/pheatmap_traj_${cor}

## 拟时序
${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--geneset_file ${geneset_file} \
--geneset_type ${geneset_type} \
--out_path ${output_dir}/celltype_plot/pheatmap_traj_${cor}

## 拟时序的顺序用自定义的
${Rscript_archr} ${scripts_dir}/run_trajectory.motifSet.orderuse.R \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--gene_order_file ${output_dir}/celltype_plot/pheatmap_traj_${cor}/gene_oder.list \
--geneset_file ${geneset_file} \
--geneset_type ${geneset_type} \
--out_path ${output_dir}/celltype_plot/pheatmap_traj_${cor}


################################################################################################
## 在germ中鉴定所有tf其对于的靶基因
################################################################################################
type=germ
${Rscript_archr} ${scripts_dir}/run_TFregulators.allTF.R \
--cor_file ${output_dir}/celltype_plot/mfuzz/cor.motif_atac-rna.tsv \
--comine_data_all_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/tf_regulators_all_${type}

## 对160个调控基因，重新做通路富集
## 每个基因的靶基因两个标准，第一个标准：cor > 0.25,LS > 80% LS;,第二个标准，新加motif enrich > 80%
type=germ
${Rscript_archr} ${scripts_dir}/run_TFregulators.allTF.pathway.R \
--gene_list_file ${output_dir}/celltype_plot/trajectory/positive/pct_0.25.list \
--input_path ${output_dir}/celltype_plot/tf_regulators_all_${type} \
--comine_data_all_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway

## 所有769个TF合并到一张表格
# 头行 
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}/*LS.tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0}' \
> ${output_dir}/celltype_plot/tf_regulators_all_${type}/germ_TF-TargetGene.tsv
# 内容
for file in `ls ${output_dir}/celltype_plot/tf_regulators_all_${type}/ | grep LS.tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}/${file} | grep -v Correlation | \
awk -F'\t' '{OFS="\t"}{print TF,$0}' TF=${tf} \
>> ${output_dir}/celltype_plot/tf_regulators_all_${type}/germ_TF-TargetGene.tsv
done
# 提取存在调控的TF-gene
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}/germ_TF-TargetGene.tsv | grep -E "YES|Correlation" \
> ${output_dir}/celltype_plot/tf_regulators_all_${type}/germ_TF-TargetGene.onlyTargetGene.tsv


## 所有160个TF合并到一张表格
# 头行 
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/*LS.tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0}' \
> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.tsv
# 内容
for file in `ls ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/ | grep LS.tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/${file} | grep -v Correlation | \
awk -F'\t' '{OFS="\t"}{print TF,$0}' TF=${tf} \
>> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.tsv
done
# 提取存在调控的TF-gene
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.tsv | grep -E "YES|Correlation" \
> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.onlyTargetGene.tsv

## 所有TF调控的基因数量合并到一张表格
echo -e "TF\tputative_targets_cor_ls_num\tputative_targets_cor_ls_enrich_num" \
> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.num.tsv
for tf in `cat ${output_dir}/celltype_plot/trajectory/positive/pct_0.25.list`
do
cor_ls_num=`cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.onlyTargetGene.tsv | \
awk -F'\t' '{if($1==tf && $11=="YES" )print}' tf=${tf} | wc -l`
cor_ls_menrich_num=`cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.onlyTargetGene.tsv | \
awk -F'\t' '{if($1==tf && $12=="YES" )print}' tf=${tf} | wc -l`
echo -e "${tf}\t${cor_ls_num}\t${cor_ls_menrich_num}" \
>> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.num.tsv
done

## 所有通路合并到一张表格,两个标准的
# 头行 
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/*ls.tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0}' \
> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.GO.tsv
# 内容
for file in `ls ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/ | grep ls.tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/${file} | grep -v gene_in_pathway | \
awk -F'\t' '{OFS="\t"}{print TF,$0}' TF=${tf} \
>> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.GO.tsv
done

## 所有通路合并到一张表格,三个标准的
# 头行 
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/*enrich.tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0}' \
> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.GO.MotifEnrich.tsv
# 内容
for file in `ls ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/ | grep enrich.tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/${file} | grep -v gene_in_pathway | \
awk -F'\t' '{OFS="\t"}{print TF,$0}' TF=${tf} \
>> ${output_dir}/celltype_plot/tf_regulators_all_${type}_pathway/germ_TF-TargetGene.GO.MotifEnrich.tsv
done

## 基于两个标准的展示每个TF的通路的热图
for geneset_file in `find ${output_dir}/celltype_plot/trajectory/positive/ | grep pct_0.25.cor_ | grep plotTrajectoryHeatmap.GEM_MM.tsv `
do
echo ${geneset_file}
type=`echo ${geneset_file} | awk -F'/' '{print $(NF-1)}'`
${Rscript_cpbd} ${scripts_dir}/tf_pathway.plot.R \
--gene_list_file ${geneset_file} \
--tf_pathway_file ${config_path}/pathway_combine.csv \
--pathway_combine_file ${output_dir}/celltype_plot/tf_regulators_all_germ_pathway/germ_TF-TargetGene.GO.tsv \
--type ${type} \
--out_path ${output_dir}/celltype_plot/tf_regulators_all_germ_pathway_plot
done

##################################################
## 计算160个TF，其motif所在peak和其它motif是否存在富集
type=germ
geneset_file=${output_dir}/celltype_plot/trajectory/positive/pct_0.25.list
${Rscript_archr} ${scripts_dir}/tf_motif.enrich.R \
--geneset_file ${geneset_file} \
--comine_data_file ${output_dir}/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/celltype_plot/trajectory/positive/pct_0.25/enrichment

## 计算160个TF，其motif是否存在协同
type=germ
cpu=20
geneset_file=${output_dir}/celltype_plot/trajectory/positive/pct_0.25.list
${Rscript_archr} ${scripts_dir}/get_motif-motif_cor.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--use_tf_file ${geneset_file} \
--cpu ${cpu} \
--cluster ${type} \
--out_path ${output_dir}/celltype_plot/trajectory/positive/pct_0.25/motif_motif
## 可视化
${Rscript_archr} ${scripts_dir}/get_motif-motif_cor.plot.R \
--cor_data_file ${output_dir}/celltype_plot/trajectory/positive/pct_0.25/motif_motif/${type}.motif_cor-chromvar.tsv \
--synergy_data_file ${output_dir}/celltype_plot/trajectory/positive/pct_0.25/motif_motif/${type}.motif_synergy-chromvar.tsv \
--out_path ${output_dir}/celltype_plot/trajectory/positive/pct_0.25/motif_motif


################################################################################################
## 分每个细胞亚群,分别鉴定motif以及peak2gene以及TF-gene的对应关系
################################################################################################
##########################################
## 鉴定motif以及peak2gene
## peak2gene默认是250kb内
## 生殖细胞
<<EOF
export type=germ
cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} | awk -F, '{print $1}' | xargs -P 5 -i sh -c '
sh ${scripts_dir}/chromvar_p2gene.sh {} ${config_path} ${type}
'
EOF

## 体细胞
## 只用peak2gene的peak
export type=somatic
cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} | awk -F, '{print $1}' | xargs -P 10 -i sh -c '
sh ${scripts_dir}/chromvar_p2gene.sh {} ${config_path} ${type}
'

## 生殖细胞
## 只用peak2gene的peak
export type=germ
cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} | awk -F, '{print $1}' | xargs -P 10 -i sh -c '
sh ${scripts_dir}/chromvar_p2gene.sh {} ${config_path} ${type}
'


## NKT细胞数量太少无法peak-gene,NKT细胞不做那一步
sh ${scripts_dir}/chromvar_p2gene.NKT.sh cluster16 ${config_path} somatic

<<EOF
## 同一细胞类型的peak在不同样本间的差异，检查是否存在差异，用于展示质控
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`

${Rscript_archr} ${scripts_dir}/peak_qc_plot.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/cell_qc/peakMA
done

## 每个细胞类型分别输出motif
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
${Rscript_archr} ${scripts_dir}/chromvar_getmotif.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--out_path ${output_dir}/motif_position
done
EOF

## 得到每个细胞的激活的peak以及motif的数量
## peak的counts>0，motif的deviations>0
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
${Rscript_archr} ${scripts_dir}/get_peak-motif_bycell.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--out_path ${output_dir}/motif_position
done

## 合并成一个文件
rm -rf ${output_dir}/motif_position/All.motif_peak.tsv
cat ${output_dir}/motif_position/cluster0.motif_peak.tsv | head -1 > ${output_dir}/motif_position/All.motif_peak.tsv
cat ${output_dir}/motif_position/*tsv | grep -v cell_type >> ${output_dir}/motif_position/All.motif_peak.tsv

## 可视化活性peak和motif的数量
${Rscript_archr} ${scripts_dir}/get_peak-motif_bycell.plot.R \
--input_file ${output_dir}/motif_position/All.motif_peak.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/motif_position

## 得到每类细胞中，基因的开放和表达的相关性以及这个基因连接的peak数量
## 胚系
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
${Rscript_archr} ${scripts_dir}/get_exp-score_peaklink.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--cluster ${cluster} \
--mean_expr_file ${output_dir}/qc_atac_v3/grem/GeneExpression.MeanByCellType.tsv \
--out_path ${output_dir}/celltype_plot/exp_score
done
## 体细胞
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
${Rscript_archr} ${scripts_dir}/get_exp-score_peaklink.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--cluster ${cluster} \
--mean_expr_file ${output_dir}/qc_atac_v3/somatic/GeneExpression.MeanByCellType.tsv \
--out_path ${output_dir}/celltype_plot/exp_score
done

## 得到每类细胞中，motif和motif间的相关性以及协同作用
## 这一步骤非常费时间
## 用已报道的motif
cat ${config_path}/Dynamic_TF.csv ${config_path}/Human_reported_TF_synergy.csv \
> ${output_dir}/celltype_plot/motif_motif/use_tf.list
cpu=10
for line in `cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell}`
do
cluster=`echo ${line} | awk -F, '{print $1}'`
${Rscript_archr} ${scripts_dir}/get_motif-motif_cor.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--use_tf_file ${output_dir}/celltype_plot/motif_motif/use_tf.list \
--cpu ${cpu} \
--cluster ${cluster} \
--out_path ${output_dir}/celltype_plot/motif_motif
done

type=germ
${Rscript_archr} ${scripts_dir}/get_motif-motif_cor.R \
--comine_data_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--use_tf_file ${output_dir}/celltype_plot/motif_motif/use_tf.list \
--cpu ${cpu} \
--cluster ${type} \
--out_path ${output_dir}/celltype_plot/motif_motif


<<EOF
最后未用
##########################################
## 用新的参数鉴定TF-gene的可能关系
## 在生殖细胞里面跑
export cor=0.5
export maxDelta=0.75
export out_path=${output_dir}/tf_regulators_${cor}_${maxDelta}

## 生殖细胞
export type=germ
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
sh ${scripts_dir}/run_TFregulators.use.sh ${cluster} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
done

## 体细胞
export type=somatic
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
sh ${scripts_dir}/run_TFregulators.use.sh ${cluster} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
done

## 链接germ和somotic到all
mkdir -p ${out_path}/all
ln -snf ${out_path}/germ/* ${out_path}/all
ln -snf ${out_path}/somatic/* ${out_path}/all


## 所有细胞类型中阳性的TF合并到一张表格
## 头行
export out_path=${output_dir}/tf_regulators_${cor}_${maxDelta}/all

cat ${out_path}/cluster0/TFregulatorPlots_All/*tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${out_path}/Positve_TF-Gene.tsv
## 具体数据行
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`
echo ${cell_type}
for file in `ls ${out_path}/${cluster}/TFregulatorPlots_All/ | grep tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${out_path}/${cluster}/TFregulatorPlots_All/${file} | grep -v putative_targets | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/Positve_TF-Gene.tsv
done
done

## 提取存在调控的TF-gene
cat ${out_path}/Positve_TF-Gene.tsv | grep -E "YES|putative_targets" \
> ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv

## 提取存在调控的TF-gene的基因列表
mkdir -p ${out_path}/gene_list
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cat ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv | grep -w ${cluster}| \
awk -F'\t' '{OFS="\t"}{print $1,$2}' | tr '\t' '\n' | sort -u \
> ${out_path}/gene_list/${cluster}.gene_interest.list
done

## 计算存在调控的TF-gene基因之间表达的相关性，用于确定TF-gene的共表达
## 内存容易溢出
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
if [ ! -f ${out_path}/expression_correlation/${cluster}_expressionCorrelation.tsv ]
then
sh ${scripts_dir}/gene_correlation.sh ${cluster} ${config_path} \
${out_path}/gene_list/${cluster}.gene_interest.list \
${out_path}/expression_correlation
fi
done

## 注释存在调控的TF-gene表达相关性
export cpu=10
sh ${scripts_dir}/tf_gene.AnnoationExpressionCor.sh ${config_path} ${cpu} ${out_path}

EOF

##########################################
## maxdelta用所有的
## 用新的参数鉴定TF-gene的可能关系
## 在生殖细胞里面跑
export cor=0.5
export maxDelta=0.75
export out_path=${output_dir}/tf_regulators_${cor}_${maxDelta}_raw

<<EOF
## 生殖细胞
export type=germ
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
sh ${scripts_dir}/run_TFregulators.use.raw.sh ${cluster} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
done
EOF

## 体细胞
<<EOF
export type=somatic
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
sh ${scripts_dir}/run_TFregulators.use.raw.sh ${cluster} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
done
EOF

## 体细胞的，只看阳性的TF
export type=somatic
cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} | awk -F, '{print $1}' | xargs -P 10 -i sh -c '
sh ${scripts_dir}/run_TFregulators.use.raw.sh {} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
'

## 生殖细胞的，所有存在表达和motif活性的都算
export type=germ
cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} | awk -F, '{print $1}' | xargs -P 10 -i sh -c '
sh ${scripts_dir}/run_TFregulators.use.allmotif.sh {} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
'


#########################
## 链接somotic到all
mkdir -p ${out_path}/all
#ln -snf ${out_path}/germ/* ${out_path}/all
ln -snf ${out_path}/somatic/* ${out_path}/all

## 所有细胞类型中阳性的TF合并到一张表格
## 头行
export out_path=${out_path}/all

cat ${out_path}/cluster0/TFregulatorPlots_All/*tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${out_path}/Positve_TF-Gene.tsv
## 具体数据行
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`
echo ${cell_type}
for file in `ls ${out_path}/${cluster}/TFregulatorPlots_All/ | grep tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${out_path}/${cluster}/TFregulatorPlots_All/${file} | grep -v putative_targets | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/Positve_TF-Gene.tsv
done
done

## 提取存在调控的TF-gene
cat ${out_path}/Positve_TF-Gene.tsv | grep -E "YES|putative_targets" \
> ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv

## 提取存在调控的TF-gene的基因列表
mkdir -p ${out_path}/gene_list
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cat ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv | grep -w ${cluster}| \
awk -F'\t' '{OFS="\t"}{print $1,$2}' | tr '\t' '\n' | sort -u \
> ${out_path}/gene_list/${cluster}.gene_interest.list
done

#########################
## 对于germ
## 所有细胞类型中阳性的TF合并到一张表格
## 头行
export out_path=${out_path}/germ

cat ${out_path}/cluster5/TFregulatorPlots_All/*tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${out_path}/Positve_TF-Gene.tsv
## 具体数据行
for line in ` cat ${config_path}/cluster_celltype.csv | grep ${germ_cell} | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`
echo ${cell_type}
for file in `ls ${out_path}/${cluster}/TFregulatorPlots_All/ | grep tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${out_path}/${cluster}/TFregulatorPlots_All/${file} | grep -v putative_targets | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/Positve_TF-Gene.tsv
done
done

## 提取存在调控的TF-gene
cat ${out_path}/Positve_TF-Gene.tsv | grep -E "YES|putative_targets" \
> ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv


########################################
## 对于之前鉴定的TF，重新跑靶基因的图和通路富集
export type=somatic
mkdir -p ${out_path}/somatic_pathway/
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`

## 提取阳性的TF列表
cat ${out_path}/somatic/${cluster}/corGIM_MM_posTFregulators.tsv | awk -F'\t' '{if($16=="YES")print $1}' | sort -u \
> ${out_path}/somatic_pathway/${cluster}.postivetf.list

${Rscript_archr} ${scripts_dir}/run_TFregulators.allTF.pathway.R \
--gene_list_file ${out_path}/somatic_pathway/${cluster}.postivetf.list \
--input_path ${out_path}/somatic/${cluster}/TFregulatorPlots_All \
--comine_data_all_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${out_path}/somatic_pathway/${cluster}
done

## 所有TF合并到一张表格
# 头行 
cat ${out_path}/somatic_pathway/cluster0/*LS.tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${out_path}/somatic_pathway/somtic_TF-TargetGene.tsv
# 内容
for file in `find ${out_path}/somatic_pathway/cluster* | grep LS.tsv`
do
echo ${file}
tf=`echo ${file} | awk -F'/' '{print $NF}' | awk -F'_' '{print $1}'`
cluster=`echo ${file} | awk -F'/' '{print $(NF-1)}'`
cell_type=`cat ${config_path}/cluster_celltype.csv | grep -w ${cluster} | awk -F, '{print $2}' | tr ' ' '_' | tr '&' '-'`
cat ${file} | grep -v Correlation | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/somatic_pathway/somtic_TF-TargetGene.tsv
done
# 提取存在调控的TF-gene
cat ${out_path}/somatic_pathway/somtic_TF-TargetGene.tsv | grep -E "YES|Correlation" \
> ${out_path}/somatic_pathway/somtic_TF-TargetGene.onlyTargetGene.tsv

#### 头行
## 所有TF调控的基因数量合并到一张表格
echo -e "cell_type\tcluster\tTF\tputative_targets_cor_ls_num\tputative_targets_cor_ls_enrich_num" \
> ${out_path}/somatic_pathway/somtic_TF-TargetGene.num.tsv
## 所有通路合并到一张表格,两个标准的
cat ${out_path}/somatic_pathway/cluster0/*ls.tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "cell_type","cluster","TF",$0}' \
> ${out_path}/somatic_pathway/somatic_TF-TargetGene.GO.tsv
## 所有通路合并到一张表格,三个标准的
cat ${out_path}/somatic_pathway/cluster0/*enrich.tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "cell_type","cluster","TF",$0}' \
> ${out_path}/somatic_pathway/somatic_TF-TargetGene.GO.MotifEnrich.tsv

#### 按照细胞类型合并
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} `
do
echo ${line}
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`

## 所有TF调控的基因数量合并到一张表格
for tf in `cat ${out_path}/somatic_pathway/${cluster}.postivetf.list`
do
cor_ls_num=`cat ${out_path}/somatic_pathway/somtic_TF-TargetGene.onlyTargetGene.tsv | grep -w ${cluster} | \
awk -F'\t' '{if($1==tf && $11=="YES" )print}' tf=${tf} | wc -l`
cor_ls_menrich_num=`cat ${out_path}/somatic_pathway/somtic_TF-TargetGene.onlyTargetGene.tsv | grep -w ${cluster} | \
awk -F'\t' '{if($1==tf && $12=="YES" )print}' tf=${tf} | wc -l`
echo -e "${cell_type}\t${cluster}\t${tf}\t${cor_ls_num}\t${cor_ls_menrich_num}" \
>> ${out_path}/somatic_pathway/somtic_TF-TargetGene.num.tsv
done

## 所有通路合并到一张表格,两个标准的
for file in `ls ${out_path}/somatic_pathway/${cluster} | grep ls.tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${out_path}/somatic_pathway/${cluster}/${file} | grep -v gene_in_pathway | \
awk -F'\t' '{OFS="\t"}{print cell_type,cluster,TF,$0}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/somatic_pathway/somatic_TF-TargetGene.GO.tsv
done

## 所有通路合并到一张表格,三个标准的
for file in `ls ${out_path}/somatic_pathway/${cluster} | grep enrich.tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${out_path}/somatic_pathway/${cluster}/${file} | grep -v gene_in_pathway | \
awk -F'\t' '{OFS="\t"}{print cell_type,cluster,TF,$0}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/somatic_pathway/somatic_TF-TargetGene.GO.MotifEnrich.tsv
done

done

################################################################################################
#### 对于sperm里面，提取存在peak2gene的peak
${Rscript_archr} ${scripts_dir}/subcell_getpeak2gene.R \
--comine_data_file ${output_dir}/subcell/cluster2/cluster2.combineRNA.motif_peak2gene.Rdata \
--out_path ${output_dir}/celltype_plot/sperm_enhancer 

## 提取bed
cat ${output_dir}/celltype_plot/sperm_enhancer/peakToGeneHeatmap_LabelClust_k25.tsv | \
awk -F'\t' '{OFS="\t"}{print $10,$11,$12}' | sed '1d' | sed 's/"//g' |\
sort -u \
> ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.bed

## 对于已报道的sperm的bed转hg38版本
cd ${output_dir}/celltype_plot/sperm_enhancer/
~/tools/StandTools/CrossMap.py bed \
~/ref/liftOver_ref/hg19ToHg38.over.chain.gz \
${output_dir}/celltype_plot/sperm_enhancer/Sperm.bed.txt \
${output_dir}/celltype_plot/sperm_enhancer/Sperm.bed.hg38.bed

#### 以下会拆分bed
## 两个bed文件取交集
${bedtools} intersect -a ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.bed \
-b ${output_dir}/celltype_plot/sperm_enhancer/Sperm.bed.hg38.bed |\
awk -F'\t' '{OFS="\t"}{print $1,$2,$3,$3-$2}' > ${output_dir}/celltype_plot/sperm_enhancer/overlap.bed

## peak2gene独立存在
${bedtools} subtract -a ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.bed \
-b ${output_dir}/celltype_plot/sperm_enhancer/Sperm.bed.hg38.bed |\
awk -F'\t' '{OFS="\t"}{print $1,$2,$3,$3-$2}' >  ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.uniq.bed

## sperm独立存在
${bedtools} subtract -a ${output_dir}/celltype_plot/sperm_enhancer/Sperm.bed.hg38.bed \
-b ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.bed |\
awk -F'\t' '{OFS="\t"}{print $1,$2,$3,$3-$2}' > ${output_dir}/celltype_plot/sperm_enhancer/Sperm_public.uniq.bed

#### 看peak的数量，不拆分bed
# -c 选项会计算 file1.bed 中每个区域与 file2.bed 中重叠区域的数量。
${bedtools} intersect -a ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.bed \
-b ${output_dir}/celltype_plot/sperm_enhancer/Sperm.bed.hg38.bed -c \
> ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.annotaionOverlap.bed

${bedtools} intersect -a ${output_dir}/celltype_plot/sperm_enhancer/Sperm.bed.hg38.bed \
-b ${output_dir}/celltype_plot/sperm_enhancer/peak2gene_sperm.bed -c \
> ${output_dir}/celltype_plot/sperm_enhancer/Sperm_public.hg38annotaionOverlap.bed

#### sperm的细胞重新聚类,共两类，并输出各自motif的值
${Rscript_archr} ${scripts_dir}/sperm_recluster.atac.R \
--comine_data_file ${output_dir}/subcell/cluster2/cluster2.combineRNA.motif_peak2gene.Rdata \
--motif_all_file ${output_dir}/celltype_plot/peak2gene/germ_mfuzz/Motif.Peak-Gene.positive.rds \
--out_path ${output_dir}/celltype_plot/sperm_recluster

#### sperm的细胞重新聚类,基于RNA,并标注atac的活性
${Rscript_archr} ${scripts_dir}/sperm_recluster.rna.R \
--cell_cluster_file ${output_dir}/celltype_plot/sperm_recluster/sperm_atac_cluster.tsv \
--comine_data_file ${output_dir}/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata \
--motif_all_file ${output_dir}/celltype_plot/peak2gene/germ_mfuzz/Motif.Peak-Gene.positive.rds \
--out_path ${output_dir}/celltype_plot/sperm_recluster

## 基于atac对sperm的重聚类结果在所有的germ细胞中标记位置
${Rscript_archr} ${scripts_dir}/sperm_recluster_show.R \
--rna_file ${output_dir}/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--cell_cluster_file ${output_dir}/celltype_plot/sperm_recluster/sperm_atac_cluster.tsv \
--out_path ${output_dir}/celltype_plot/sperm_recluster

## 基于atac对sperm的重聚类结果分别计算peak2gene
${Rscript_archr} ${scripts_dir}/sperm_recluster.peak2gene.R \
--comine_data_file ${output_dir}/subcell/cluster2/cluster2.combineRNA.motif_peak2gene.Rdata \
--cell_cluster_file ${output_dir}/celltype_plot/sperm_recluster/sperm_atac_cluster.tsv \
--out_path ${output_dir}/celltype_plot/sperm_recluster

################################################################################################
## 对所有细胞都尝试重新聚类,基于ATAC
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-'`
do
echo ${line}
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`

${Rscript_archr} ${scripts_dir}/sperm_recluster.atac.all.R \
--cluster ${cluster} \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--motif_all_file ${output_dir}/celltype_plot/peak2gene/germ_mfuzz/Motif.Peak-Gene.positive.rds \
--out_path ${output_dir}/celltype_plot/recluster_allcell &
done

## 基于atac对所有的重聚类结果在所有的细胞中标记位置,germ
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} `
do
echo ${line}
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`

${Rscript_archr} ${scripts_dir}/sperm_recluster_show_all.R \
--cluster ${cluster} \
--rna_file ${output_dir}/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--cell_cluster_file ${output_dir}/celltype_plot/recluster_allcell/${cluster}_atac_cluster.tsv \
--out_path ${output_dir}/celltype_plot/recluster_allcell & 
done
# somatic
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} `
do
echo ${line}
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`

${Rscript_archr} ${scripts_dir}/sperm_recluster_show_all.R \
--cluster ${cluster} \
--rna_file ${output_dir}/qc_atac_v3/somatic/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac_v3/somatic/testis_combined_peak.combineRNA.qc.Rdata \
--cell_cluster_file ${output_dir}/celltype_plot/recluster_allcell/${cluster}_atac_cluster.tsv \
--out_path ${output_dir}/celltype_plot/recluster_allcell & 
done


################################################################################################
## 所有已报道的关键的TF在motif水平、其本身开放程度以及表达水平进行染色
for gene in `cat ${config_path}/Human_reported_TF2.csv | grep -v gene | awk -F, '{print $1}' | tr '/' '\n'`
do
${Rscript_archr} ${scripts_dir}/plotEmbedding_TFregulators.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--comine_data_all_file ${output_dir}/qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--gene ${gene} \
--out_path ${output_dir}/report_tf
done

## 下面快,基因读进去
${Rscript_archr} ${scripts_dir}/plotEmbedding_TFregulators.readlist.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--comine_data_all_file ${output_dir}/qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--gene_list_file ${config_path}/Human_reported_TF2.csv \
--out_path ${output_dir}/report_tf

## 只看germ
${Rscript_archr} ${scripts_dir}/plotEmbedding_TFregulators.readlist.R \
--rna_file ${output_dir}/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata \
--comine_data_all_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--gene_list_file ${config_path}/Human_reported_TF2.csv \
--out_path ${output_dir}/report_tf_germ
${Rscript_archr} ${scripts_dir}/plotEmbedding_TFregulators.onlyRNA.readlist.R \
--rna_file ${output_dir}/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata \
--gene_list_file ${config_path}/Human_reported_TF2.csv \
--out_path ${output_dir}/report_tf_germ

## 只看rna
gene_list=(PIWIL4 GFRA1 L1TD1 UTF1 EGR4 ETV5)
for gene in ${gene_list[@]}
do
${Rscript_archr} ${scripts_dir}/plotEmbedding_TFregulators.onlyRNA.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--gene ${gene} \
--out_path ${output_dir}/report_tf
done


## 所有已报道的关键的TF其在germ中的调控情况
for gene in `cat ${config_path}/Human_reported_TF2.csv | grep -v gene | awk -F, '{print $1}' | tr '/' '\n'`
do
${Rscript_archr} ${scripts_dir}/plot_peak-gene.R \
--comine_data_all_file ${output_dir}/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--gene ${gene} \
--out_path ${output_dir}/report_tf
done

## 所有已报道的转录因子的聚类,等画
for gene in `cat ${config_path}/dureplicate_reported_TF.csv | grep -v gene | awk -F, '{print $1}'`
do
${Rscript_archr} ${scripts_dir}/plot_peak-gene.R \
--comine_data_all_file ${output_dir}/qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--gene ${gene} \
--out_path ${output_dir}/report_tf
done

## 所有marker基因关键的TF其在germ和somatic中的调控情况
type_list=(germ somatic)
for type in ${type_list[@]}
do
for gene in `cat ${config_path}/marker_gene_peak.csv | grep -v Symbol | awk -F, '{print $2}'`
do
${Rscript_archr} ${scripts_dir}/plot_peak-gene.R \
--comine_data_all_file ${output_dir}/qc_atac_v3/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--gene ${gene} \
--out_path ${output_dir}/celltype_plot/marker_gene_peak_${type}
done
done

################################################################################################
#### 标记并展示已知的重要基因在rna及atac的热图
################################################################################################
## 基因分为总的，TF以及nonTF
## 展示按照基因的表达聚类以及按照基因的固定顺序
${Rscript_archr} ${scripts_dir}/pheatmap_knowngene.R \
--rna_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--comine_data_all_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--known_gene_file ${config_path}/SCOS_known_genes.csv \
--out_path ${output_dir}/celltype_plot/pheatmap_knowngene



################################################################################################
#### 细胞通讯
################################################################################################
## 按照两种分类，生殖细胞9类和所有体细胞做
## 另外一类是
${Rscript_cpbd} ${scripts_dir}/cellphonedb.R \
--rna_file ${output_dir}/qc_atac_v3/all/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/celltype_plot/cpdb

## python运行
mkdir -p ${output_dir}/celltype_plot/cpdb/divide
mkdir -p ${output_dir}/celltype_plot/cpdb/combine

## 生殖细胞为9类细胞
${python_cpbd} ${scripts_dir}/cellphonedb.divide.py
## 生殖细胞为3类细胞
${python_cpbd} ${scripts_dir}/cellphonedb.combine.py

## 可视化
type=divide
${Rscript_cpbd} ${scripts_dir}/cellphonedb_ktplots.R \
--pvals_file ${output_dir}/celltype_plot/cpdb/${type}/statistical_analysis_pvalues_germ_somatic.txt \
--means_file ${output_dir}/celltype_plot/cpdb/${type}/statistical_analysis_significant_means_germ_somatic.txt \
--gene_list_file ${output_dir}/celltype_plot/cpdb/${type}/show_gene.list \
--type ${type} \
--out_path ${output_dir}/celltype_plot/cpdb/${type}




################################################################################################
## 共享数据
mkdir -p ${output_dir}/shareData/peak-gene
mkdir -p ${output_dir}/shareData/motif/cisbp
mkdir -p ${output_dir}/shareData/motif/jaspar

cd ${output_dir}/shareData

## 细胞类型注释
cp -rf ${config_path}/cluster_celltype.csv ${output_dir}/shareData/

## 转录组数据
ln -snf ${input_dir}/testis_combined.annotationCellType.Rdata ${output_dir}/shareData/

## atac包含peak的信息
ln -snf ${output_dir}/atac_res/testis_combined_peak.Rdata ${output_dir}/shareData/

## peak-gene
for file in `find ${output_dir}/cluster_all_result | grep correlation_magic.Rdata`
do
ln -snf ${file} ${output_dir}/shareData/peak-gene
done

## motif
for file in `find ${output_dir}/motif | grep rda | grep cisbp`
do
ln -snf ${file} ${output_dir}/shareData/motif/cisbp
done
for file in `find ${output_dir}/motif | grep rda | grep jaspar`
do
ln -snf ${file} ${output_dir}/shareData/motif/jaspar
done